expit <- function(x) exp(x)/(1+exp(x))

#############################
#############################
#############################

sparsemix_media <- function(votes.mat,missing.mat,max.loop=500,
                      burnin=50, thin=1,empir.cdf=NULL,cutoff.seq=0,
                      save.curr="UDV_curr",
                      save.each=FALSE,thin.save=25,max.optim=300, verbose=FALSE, seed){

    if (!missing(seed)){ 
        set.seed(seed)
    }
    ##votes.mat=dtm2;missing.mat=dtm2*0;max.loop=500; burnin=50; thin=1;empir.cdf=NULL;cutoff.seq=seq.vec;save.curr="UDV_curr";save.each=FALSE;thin.save=25;max.optim=300

    n<-dim(votes.mat)[1]
    k<-dim(votes.mat)[2]

    votes.mat<-as.matrix(votes.mat)

    row.type<-apply(votes.mat,1,FUN=function(x)  (max(x)<=1))
    row.type[row.type==TRUE]<-"bin"
    row.type[row.type==FALSE]<-"count"
    if(verbose){
        print(table(row.type))
    }
    sigma<-1

    ##Initialize D, lambda, Theta, Z
    votes.mat.init<-t(apply(votes.mat,1,FUN=function(x) x/sd(x)))
    votes.mat.init[!is.finite(votes.mat.init)]<-0
    mu.r<-rowMeans(sigma*votes.mat.init)
    mu.c<-colMeans(sigma*votes.mat.init)
    mu.grand<- mean(sigma*votes.mat.init)

    ones.r<-rep(1,k)
    ones.c<-rep(1,n)


    svd.mat<- sigma*votes.mat.init #- ones.c%*%t(mu.c)-mu.r%*%t(ones.r)+mu.grand
    svd.mat<- svd.mat-ones.c%*%t(colMeans(svd.mat))-rowMeans(svd.mat)%*%t(ones.r)+mean(svd.mat)


    Dtau<-rep(1, dim(votes.mat)[1])
    Theta.last<-Z.next<-Z.last<-0*votes.mat
    Theta.last<-Z.last<-U.last<-0*votes.mat

    ##Containers
    U.run<-U.run.2<-matrix(NA, nr=dim(votes.mat)[1], nc=max.loop+burnin)
    V.run<-V.run.2<-matrix(NA, nr=dim(votes.mat)[2], nc=max.loop+burnin)
    D.run<-D.mean.run<-matrix(NA, nr=min(n,k), nc=max.loop+burnin)
    accept.run<-NULL
    lambda.run<-NULL
    dev.run<-NULL
    V.average<-U.average<-0
    Z.run<-Z.mode.run<-0
    if(sum(row.type=="bin")>2){
	Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==1]<-dnorm(0)/pnorm(.5)
	Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==0]<- -dnorm(0)/pnorm(.5)
	Theta.last[row.type=="bin",][missing.mat[row.type=="bin",]==1]<- 0}

    dubcent<-function(X){
        colmean.vec<-colMeans(X,na.rm=T)
        rowmean.vec<-rowMeans(X,na.rm=T)
        for(i in 1:dim(X)[1]) for (j in 1:dim(X)[2])
            X[i,j]<-X[i,j]-rowmean.vec[i]-colmean.vec[j]
        X-mean(X)
    }

    if(sum(row.type=="count")>2){
	Theta.last[row.type=="count",]<- (log(1*(votes.mat[row.type=="count",]==0)+votes.mat[row.type=="count",]))
        ##Theta.last[row.type=="count",]<-(Theta.last[row.type=="count",]-mean(Theta.last[row.type=="count",]))/sd(Theta.last[row.type=="count",])
	Theta.last[row.type=="count",][missing.mat[row.type=="count",]==1]<- 0
        ##Theta.last[row.type=="count",]<- 0
    }

    svd0<-svd(svd.mat)
    lambda.shrink<-lambda.lasso<-rep(median((2/(svd0$d))^.5),min(n,k))[1]

    params.init<-c(0,0)#log(c(mean(votes.mat)))
    proposal.sd<-0.1
    step.size<-scale.sd<-1

#############################
#############################
#############################
###########Start loop here
    for(i.gibbs in 1:(max.loop+burnin)){

        ## Step 1: Make Z-star
        Z.call<-make.Z_media(Theta.last.0=Theta.last, votes.mat.0=votes.mat,
                       row.type.0=row.type,  n0=n ,k0=k, iter.curr=i.gibbs,
                       empir=empir.cdf,cutoff.seq.0=cutoff.seq,missing.mat.0=missing.mat,
                       lambda.lasso=lambda.lasso,params=params.init,
                       proposal.sd=proposal.sd,scale.sd=scale.sd,max.optim=max.optim,
                       step.size=step.size)

        Z.next<-Z.call$Z
        params.init<-Z.call$params
        accept.run[i.gibbs]<-Z.call$accept
        dev.run[i.gibbs]<-Z.call$prob
        params.init<-Z.call$param
        proposal.sd<-Z.call$proposal.sd
        step.size<-Z.call$step.size
        
        ## Step 2: Update UDV and components therein
        UDV.all <- update.UDV_media(
            Z.next.0=Z.next,
            k0=k, n0=n, lambda.lasso.0=lambda.lasso,lambda.shrink.0=lambda.shrink,
            votes.mat.0=votes.mat, Dtau.0=Dtau, iter.curr=i.gibbs,
            row.type.0=row.type, missing.mat.0=missing.mat,vote.wt.0=1)
        
        last.Z<-Z.call
        last.UDV<-UDV.all

        if(i.gibbs>(max.optim+50)&i.gibbs%%50==1){
            seq.check<-(i.gibbs-1):(i.gibbs-20)
            if(mean(accept.run[(i.gibbs-1):(i.gibbs-20)])>.9 ) step.size<-step.size*1.25
            if(mean(accept.run[(i.gibbs-1):(i.gibbs-20)])<.1 ) step.size<-step.size*.8
            print("acceptance percentage")
            print(mean(accept.run[(i.gibbs-1):(i.gibbs-20)]))
        }


        ## Update inputs to UDV.all
        Dtau=UDV.all$Dtau
        lambda.lasso<-UDV.all$lambda.lasso
        lambda.shrink<-UDV.all$lambda.shrink
        Theta.last<-UDV.all$Theta.last

        if(verbose){
            ##Report output
            if(i.gibbs%%1==0){
                ##print(table(sign(Theta.last),votes.mat))
                print("range theta")
                print(range(Theta.last))
                print("iter")
                print(i.gibbs)
                ##print("rms")
                ##print(mean((Theta.last-Z.next)^2))
                ##print(c("sigma",sigma))
                print("trunc")
                print(UDV.all$D.trunc)
                print("post")
                print(UDV.all$D.post[1:length(UDV.all$D.trunc)])
                print("lambda.mean")
                print(UDV.all$lambda.lasso)
                ##print("lambda.mode")
                ##print(UDV.all$lambda.shrink)
            } ##Closes printing of results
        }
        
        ##Gather results
        lambda.run[i.gibbs]<-UDV.all$lambda.lasso[1]
        U.run[,i.gibbs]<-UDV.all$U.next[,1]
        U.run.2[,i.gibbs]<-UDV.all$U.next[,2]
        V.run[,i.gibbs]<-UDV.all$V.next[,1]
        V.run.2[,i.gibbs]<-UDV.all$V.next[,2]
        D.run[,i.gibbs]<-UDV.all$D.trunc[1:min(n,k)]
        D.mean.run[,i.gibbs]<-UDV.all$D.post[1:min(n,k)]

        if(i.gibbs>burnin & i.gibbs %% thin == 0) {
            Z.run<- Z.run + Theta.last
            Z.mode.run<- Z.mode.run + UDV.all$Theta.mode
            U.average<-U.average+UDV.all$U.next
            V.average<-V.average+UDV.all$V.next
        }
        ##if(i.gibbs %%100 ==0) print(i.gibbs)
    } #Closes out for loop
    
    ## svd.UDV<-svd(dubcent(Z.run))
    
    out <- list("D.sparse"=D.run[,-c(1:burnin)], "D.mean"=D.mean.run[,-c(1:burnin)],
                "Udim1"=U.run[,-c(1:burnin)], "Udim2"=U.run.2[,-c(1:burnin)],
                "Vdim1"=V.run[,-c(1:burnin)], "Vdim2"=V.run.2[,-c(1:burnin)],
                "lambda.lasso"=lambda.run[-c(1:burnin)], "Z" = Z.run/max.loop,
                "Z.mode"= Z.mode.run/max.loop,
                "fitted" = pnorm(Z.mode.run/max.loop), "U.all"=U.average/max.loop,
                "V.all"=V.average/max.loop,"accept"=accept.run,"dev.run"=dev.run)
                ## "fitted" = pnorm(Z.mode.run/max.loop), "U.all"=svd.UDV$u, "V.all"=svd.UDV$v,"accept"=accept.run,"dev.run"=dev.run)
    
}


################################
################################
##Three functions in here:
#### 1) test.gamma for finding the ML estimate for the grand lambda
#### 2) make.Z for converting a matrix to the Z-scale
#### 3) update.UDV for updating the ideal point estimates
#### *) Small additonal ones at the end 
################################
################################

################################
################################
##1) Finding gamma
test.gamma.pois<-function(gamma.try,Theta.last.0=Theta.last[row.type=="count",],votes.mat.0=votes.mat[row.type=="count",],
emp.cdf.0,cutoff.seq=NULL){
	
	gamma.try<-exp(gamma.try)
	votes.mat<-votes.mat.0
	taus.try<-NULL
	count.seq<-cutoff.seq
	if(length(cutoff.seq)==0) count.seq<-seq(-1,max(votes.mat)+2,1)
	taus.try<-count.seq*0

	analytic.cdf<-count.seq
	a.int<-qnorm(mean(votes.mat==0))
	taus.try[count.seq>=0]<-(a.int+gamma.try[1]*count.seq[count.seq>=0]^(gamma.try[2]))
	taus.try[count.seq<0]<- -Inf
	taus.try<-sort(taus.try)
	taus.try[1]<--Inf	
	find.pnorm<-function(x){
		a<-x[1]
		b<-x[2]
		coefs<-c( -1.82517672, 0.51283415, -0.81377290, -0.02699400, -0.49642787, -0.33379312, -0.24176661, 0.03776971)
		x<-c(1, a, b, a^2,b^2, log(abs(a-b)),log(abs(a-b))^2,a*b)
		(sum(x*coefs))
	}
	lik<-((pnorm(taus.try[votes.mat+2 ]-Theta.last.0)-pnorm(taus.try[votes.mat+1 ]-Theta.last.0)) )
	log.lik<-log(lik)
	
	which.zero<-which(lik==0)
	a0<-taus.try[votes.mat[which.zero]+2]-Theta.last.0[which.zero]
	b0<-taus.try[votes.mat[which.zero]+1]-Theta.last.0[which.zero]
	log.lik[which.zero]<-apply(cbind(a0,b0),1,find.pnorm)
	log.lik[is.infinite(lik)&lik>0]<-max(log.lik[is.finite(lik)])
	log.lik[is.infinite(lik)&lik<0]<-min(log.lik[is.finite(lik)])

	thresh<-min(-1e30, min(log.lik[is.finite(log.lik)],na.rm=TRUE))
	log.lik[log.lik<thresh]<-thresh
	dev.out<--2*sum(log.lik,na.rm=TRUE)+sum(log(gamma.try)^2)
	
	
	return(list("deviance"=dev.out,"tau"=taus.try))

}	



################################
################################
##2) Converting a matrix to the Z scale
	make.Z_media<-function(
			Theta.last.0=Theta.last, 
			votes.mat.0=votes.mat,row.type.0=row.type, 
			n0=n, k0=k, params=NULL, iter.curr=0,empir=NULL,cutoff.seq.0=NULL,missing.mat.0=NULL,lambda.lasso,proposal.sd,scale.sd, max.optim,step.size
			){
				
		Theta.last<-Theta.last.0;votes.mat<-votes.mat.0;
		row.type<-row.type.0; n<-n0; k<-k0
		cutoff.seq<-cutoff.seq.0
		sigma<-1
		row.type<-row.type.0; votes.mat<-votes.mat.0
		Z.next<-matrix(NA,nr=n,nc=k)
		missing.mat<-missing.mat.0
		i.gibbs<-iter.curr
		print(i.gibbs)
		print(iter.curr)
	

	
			if(sum(row.type=="bin")>2){
				
			Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==1]<-rtnorm(sum(votes.mat[row.type=="bin",]==1), mean=sigma^.5*Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==1], lower=0, sd=sigma^.5)
		Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==0]<-rtnorm(sum(votes.mat[row.type=="bin",]==0), mean=sigma^.5*Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==0], upper=0 , sd=sigma^.5)
		Z.next[row.type=="bin",][missing.mat[row.type=="bin",]==1]<-rnorm(sum(missing.mat[row.type=="bin",]==1), mean=sigma^.5*Theta.last[row.type=="bin",][missing.mat[row.type=="bin",]==1], sd=sigma^.5)
		
		
		pars.max<-1
		accept.out<-prob.accept<-1

				}

		if(sum(row.type=="count")>2){
		#begin type = "pois"


	num.sample.count<-sum(row.type=="count")*k
	accept.out<-prob.accept<-NA

	dev.est<-function(x) test.gamma.pois(x,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	dev.est1<-function(x) test.gamma.pois(c(x,params[2]),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	dev.est2<-function(x) test.gamma.pois(c(params[1],x),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	
	
	#range.opt<-c(params-1,params+1)
	range.opt<-rbind(params-1,params+1)
	grad.est<-function(x,del=.01){
		grad1<-(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		grad2<-(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		grad.all<-c(grad1,grad2)
		grad.all[!is.finite(grad.all)]<-0
		grad.all
		}
	grad.est.1<-function(x,del=.01){
		grad1<-(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		#grad2<-0#(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		#c(grad1,grad2)
		grad1[!is.finite(grad1)]<-0
		grad1
		}
	grad.est.2<-function(x,del=.01){
		#grad1<-0#(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		grad2<-(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		grad2[!is.finite(grad2)]<-0
		#c(grad1,grad2)
		grad2
		}

	grad.est.12<-function(x,del=.001){
		#grad1<-0#(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		#grad12<-(dev.est(x+del)-dev.est(x-del))/(8*del)
		grad12<-dev.est(x+del)-dev.est(x+del*c(1,-1))-dev.est(x+del*c(-1,1))+dev.est(x-del)
		grad12<-grad12/(8*del^2)
		grad12[!is.finite(grad12)]<-0
		grad12
		#c(grad1,grad2)
		}
		

	if(iter.curr<=max.optim){
	#gamma.opt<-optimize(dev.est,
	#	lower=range.opt[1],upper=range.opt[2])
	gamma.opt.1<-optimize(dev.est1,
			lower=range.opt[1,1],upper=range.opt[2,1])
	params[1]<-gamma.opt.1$minimum
	gamma.opt.2<-optimize(dev.est2,
		lower=range.opt[1],upper=range.opt[2])
	params[2]<-gamma.opt.2$minimum
	gamma.next<-pars.max<-params
	gamma.opt<-list(gamma.opt.1,gamma.opt.2)
	#gamma.opt<-optim(params,dev.est,method="Nelder-Mead")
		#pars.max<-exp(gamma.opt$min)
	#	pars.max<-exp(gamma.opt$par)
	#	proposal.sd<-1
	#	gamma.next<-pars.max
		#proposal.sd<-abs(params-gamma.opt$min)
		} else {
		##Hamiltonian step
		##Code adapted from Neal, Handbook of MCMC
		gamma.opt<-"Hamiltonian Step"
		num.HMC<-10
		#q.run<-p.run<-matrix(NA,nr=2,nc=(num.HMC*2))
		if(i.gibbs==1) step.size<-1
		if(i.gibbs==20) step.size<-.5
		if(i.gibbs==50) step.size<-.2
		if(i.gibbs==100) step.size<-.1	
#		print(accept.run)
		#if(i.gibbs>150& i.gibbs%%20==1){
		#	if(mean(accept.run[(i.gibbs-1):(i.gibbs-50)])>.6) step.size<-step.size*1.25
		#	if(mean(accept.run[(i.gibbs-1):(i.gibbs-50)])<.2) step.size<-step.size*.8
		#}	

		
		step.size<-min(step.size,2)
		
		print("step size")
		print(step.size)
		
		q0<-q<-params
		p0<-p<-rnorm(length(q),0,1)
		q
		p
		num.HMC<-10
		q.run<-p.run<-matrix(NA,nr=2,nc=(num.HMC*2))
		for(i.HMC in 1:(num.HMC*2)){
		if(i.HMC%%3==1){
		hess.mat<-matrix(NA,nr=2,nc=2)
		hess.mat[1,1]<-grad.est.1(q)
		hess.mat[1,2]<-hess.mat[2,1]<-grad.est.12(q)
		hess.mat[2,2]<-grad.est.2(q)
		hess.mat<-.5*(hess.mat+t(hess.mat))
		epsilon<-ginv(hess.mat)*step.size/10
		}
		q<-q+epsilon%*%p
		p<-p-epsilon%*%grad.est(q)
		}
		q<-q+epsilon%*%p
		p<-p-epsilon%*%grad.est(q)/2
		p<- -p
		
		current_U<-dev.est(q0)/2
		current_K<-sum(p0^2)/2
		proposed_U<-dev.est(q)/2
		proposed_K<-sum(p^2)/2
		proposal.dev<-proposed_U+proposed_K
		current.dev<-current_U+current_K
		gamma.cand<-q
		print("Log Probs")
		print(c(proposal.dev,current.dev))
		if(exp(current.dev-proposal.dev)>runif(1)){
			##Do if accept
			pars.max<-gamma.next<-(gamma.cand)
			print("accept")
			accept.out<-prob.accept<-1
			} else {
			##Do if reject
			pars.max<-gamma.next<-(params)
			print("reject")
			accept.out<-prob.accept<-0
			}
		}
		print("Proposal sd")
		print(proposal.sd)
		taus<-test.gamma.pois(gamma.next,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$tau
		#for(i.tau in 2:length(taus)) if(taus[i.tau]==taus[i.tau-1]) taus[i.tau]<-taus[i.tau]+.001
		#print("SumTau")
		#print(sum(taus[-1]))
		taus<-sort(taus)
		taus[1]<--Inf
		print("Range of gamma")
		print(range.opt)
		print(gamma.opt)
		print(range(taus))
		
		
		Z.next.temp<-rtruncnorm(num.sample.count, mean=Theta.last[row.type=="count",], a=taus[votes.mat[row.type=="count",]+1 ],b=taus[votes.mat[row.type=="count",]+2 ] )

		
		Z.next[row.type=="count",]<-matrix(Z.next.temp,nrow=sum(row.type=="count"))



		}
		
	return(list("Z.next"=Z.next,"params"=(pars.max),"accept"=accept.out,"prob"=prob.accept,"proposal.sd"=proposal.sd,"step.size"=step.size))
	}##Closes out make.Z function


################################
################################
##2) Converting a matrix to the Z scale
	update.UDV_media<-function(
	Z.next.0=Z.next,
	k0=k, n0=n, lambda.lasso.0=lambda.lasso,lambda.shrink.0=lambda.shrink,
	Dtau.0=Dtau,
	votes.mat.0=votes.mat, iter.curr=0,row.type.0,missing.mat.0=missing.mat,vote.wt.0=1
	){
	vote.wt<-vote.wt.0
	missing.mat<-missing.mat.0
	Dtau<-Dtau.0;
	Z.next<-Z.next.0;
	votes.mat<-votes.mat.0;
	n<-n0; k<-k0; 
	lambda.lasso<-lambda.lasso.0
	row.type <- row.type.0
	#cols.impute<-cols.impute.0

	#Declare some vectors
	sigma<-1
	ones.r<-rep(1,k0)
	ones.c<-rep(1,n0)

	#Update intercepts
	mu.r<-rowMeans(Z.next)#-ones.c%*%t(mu.c)-Theta.last.0+mu.grand)
	mu.r<-mu.r*n/(n+1)+rnorm(length(mu.r),sd=1/k)

	mu.c<-colMeans(Z.next)#-mu.r%*%t(ones.r)-Theta.last.0+mu.grand)
	mu.c<-mu.c*k/(k+1)+rnorm(length(mu.c),sd=1/n)
	mu.grand<-mean(Z.next)
	mu.grand<-mu.grand*(n*k)/(n*k+1)+rnorm(1,sd=1/(n*k))
	
	mean.mat<-ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand
	
	if(length(unique(row.type))==2){
		print("Two Means Being Used")
		mean.c.mat<-matrix(NA,nrow=n,ncol=k)
		mu.c.bin<-colMeans(Z.next[row.type=="bin",])
		mu.c.count<-colMeans(Z.next[row.type=="count",])
		mu.c.bin<-mu.c.bin*k/(k+1)+rnorm(length(mu.c.bin),sd=1/sum(row.type=="bin"))
		mu.c.count<-mu.c.count*k/(k+1)+rnorm(length(mu.c.count),sd=1/sum(row.type=="count"))
		mean.c.mat[row.type=="bin",]<-ones.c[row.type=="bin"]%*%t(mu.c.bin)
		mean.c.mat[row.type=="count",]<-ones.c[row.type=="count"]%*%t(mu.c.count)
	
		mean.grand.mat<-matrix(NA,nrow=n,ncol=k)
		mean.r.mat<-matrix(NA,nrow=n,ncol=k)

		mean.grand.mat[row.type=="bin",c(1,2)]<-mean(Z.next[row.type=="bin",c(1,2)])+rnorm(1,sd=1/(sum(row.type=="bin")*k))
		mean.grand.mat[row.type=="count",c(1,2)]<-mean(Z.next[row.type=="count",c(1,2)])+rnorm(1,sd=1/(sum(row.type=="count")*k))

		mean.grand.mat[row.type=="bin",-c(1,2)]<-mean(Z.next[row.type=="bin",-c(1,2)])+rnorm(1,sd=1/(sum(row.type=="bin")*k))
		mean.grand.mat[row.type=="count",-c(1,2)]<-mean(Z.next[row.type=="count",-c(1,2)])+rnorm(1,sd=1/(sum(row.type=="count")*k))

	mu.r.1<-rowMeans(Z.next[,c(1:3)])#-ones.c%*%t(mu.c)-Theta.last.0+mu.grand)
	mu.r.1<-mu.r.1*n/(n+1)+rnorm(length(mu.r),sd=1/k)

	mu.r.2<-rowMeans(Z.next[,-c(1:3)])#-ones.c%*%t(mu.c)-Theta.last.0+mu.grand)
	mu.r.2<-mu.r.2*n/(n+1)+rnorm(length(mu.r),sd=1/k)

	mean.r.mat[,1:3]<-mu.r.1%*%t(ones.r[1:3])
	mean.r.mat[,-c(1:3)]<-mu.r.2%*%t(ones.r[-c(1:3)])

	mean.mat<-mean.c.mat+mean.r.mat-mean.grand.mat
	
	}

	Z.starstar<-svd.mat<- Z.next- mean.mat



	

	#Take svd, give each column an sd of 1 (rather than norm of 1)
	num.zeroes<-1#colMeans(votes.mat!=0)

	
	save(svd.mat,file="svd.mat")
	
	svd.mat[is.na(svd.mat)]<-0

	svd.mat.0<-svd.mat
	
	wts.dum<-rep(1,nrow(svd.mat))
	wts.dum[row.type=="bin"]<-1/sum(1-missing.mat[row.type=="bin",])^.5*vote.wt/(1+vote.wt)
	wts.dum[row.type=="count"]<-1/sum(1-missing.mat[row.type=="count",])^.5*vote.wt/(1+vote.wt)


	wts.dum<-wts.dum/mean(wts.dum)
	
	svd.dum.0<-svd.dum<-svd(svd.mat*wts.dum)

	svd.dum$u[is.na(svd.dum$u)|is.infinite(svd.dum$u)]<-0
	svd.dum$d[is.na(svd.dum$d)|is.infinite(svd.dum$d)]<-0
	svd.dum$v[is.na(svd.dum$v)|is.infinite(svd.dum$v)]<-0


	svd.dum$d<-(t(svd.dum$u)%*%svd.mat%*%svd.dum$v)
	which.rows<-which(rowMeans(svd.dum$d^2)^.5<1e-4)
	which.cols<-which(colMeans(svd.dum$d^2)^.5<1e-4)
	svd.dum$d[which.rows,]<-rnorm(length(svd.dum$d[which.rows,]),sd=.001)
	svd.dum$d[which.cols,]<-rnorm(length(svd.dum$d[which.cols,]),sd=.001)
	svd2<-svd(svd.dum$d)
	svd2$u[is.na(svd2$u)|is.infinite(svd2$u)]<-0
	svd2$d[is.na(svd2$d)|is.infinite(svd2$d)]<-0
	svd2$v[is.na(svd2$v)|is.infinite(svd2$v)]<-0
	
	svd.dum$u<-svd.dum$u%*%(svd2$u)
	svd.dum$v<-svd.dum$v%*%(svd2$v)

	svd.dum$d<-svd2$d
	
	svd0<-svd.dum
	
	svd0$v<-t(t(svd0$u)%*%svd.mat.0)
	svd0$v<-apply(svd0$v,2,FUN=function(x) x/sum(x^2)^.5)
	svd0$d<-diag(t(svd0$u)%*%svd.mat.0%*%svd0$v)
	sort.ord<-sort(svd0$d,ind=T,dec=T)$ix
	svd0$u<-svd0$u[,sort.ord]
	svd0$v<-svd0$v[,sort.ord]
	svd0$d<-svd0$d[sort.ord]
	svd0$u<-svd0$u*(n-1)^.5
	svd0$v<-svd0$v*(k-1)^.5
	svd0$d<-svd0$d*((n-1)*(k-1))^-.5
	Theta.last.0<-svd.mat
	Theta.last<-Theta.last.0+(ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand)

	#Update d; follows from Blasso and DvD
	Y.tilde<-as.vector(svd.mat)
	if(n>length(Dtau)) Dtau[(length(Dtau)+1):n]<-1
	A<- (n*k)*diag(n)+diag(as.vector(Dtau^(-1)))
	#if(n>k) for(i.A in (k+1):n) A[i.A,i.A]<-(Dtau^-1)[i.A]*0+1e20
	#gA<-ginv(A)
	gA<-A*0+NA
	gA[1:min(n,k),1:min(n,k)]<-ginv(A[1:min(n,k),1:min(n,k)])
	gA[is.na(gA)]<-0
	XprimeY<-sapply(1:min(n,k), FUN=function(i, svd2=svd0, Z.use=svd.mat) sum(Z.use*(svd2$u[,i]%*%t(svd2$v[,i]))))
	if(length(XprimeY)<dim(gA)[1]) XprimeY[(length(XprimeY)+1) :dim(gA)[1]]<-0
	D.post.mean<- as.vector(gA%*%XprimeY)
	D.post.var.2<-gA

##Sample D and reconstruct theta, putting intercepts back in
	D.post<-rep(NA,length(D.post.mean))
	D.post[1:min(n,k)]<-as.vector(mvrnorm(1, mu=D.post.mean[1:min(n,k)], D.post.var.2[1:min(n,k),1:min(n,k)] ) )/sigma^.5
	D.post[is.na(D.post)]<-0
	abs.D.post<-abs(D.post)

##Calculate MAP and mean estimate
	D.trunc<-pmax(svd0$d-lambda.lasso.0,0)
	
	U.last<-svd0$u
	V.last<-svd0$v

	sample.UV<-TRUE
	if(iter.curr<=50) sample.UV<-FALSE
	if(sample.UV==TRUE){
	U.mean.next<-0*U.last
	V.mean.next<-0*V.last
##Construct U and V
#prior var of 2
	for(l in 1:dim(U.last)[2]){
		prior.var<-.25#/D.post[l]
	Y.star.l<- svd.mat-U.last[,-l]%*%diag(D.post[1:min(n,k)][-l])%*%t(V.last[,-l])
	U.mean.next<-rowSums(abs.D.post[l]*(Y.star.l%*%V.last[,l]))/
						(abs.D.post[l]^2*sum(V.last[,l]^2)+prior.var)
							U.var.next<-sigma^2/(abs.D.post[l]^2*sum(V.last[,l]^2) + prior.var)
	U.next.temp<-U.mean.next+rnorm(length(U.mean.next),0, sd=U.var.next^.5)
	U.next.temp<-	U.next.temp-mean(U.next.temp)
	}
	
	U.last[is.na(U.last)]<-0
	
	for(l in 1:dim(V.last)[2]){
		prior.var<-.25#/D.post[l]
	Y.star.l<- svd.mat-U.last[,-l]%*%diag(D.post[1:min(n,k)][-l])%*%t(V.last[,-l])
	V.mean.next<-colSums(abs.D.post[l]*(t(U.last[,l])%*%Y.star.l))/
						(abs.D.post[l]^2*sum(U.last[,l]^2)+prior.var)
					V.var.next<-sigma^2/
								(abs.D.post[l]^2*sum(U.last[,l]^2)+prior.var)
	V.next.temp<-V.mean.next+rnorm(length(V.mean.next),0, sd=V.var.next^.5)
	V.next.temp<-	V.next.temp-mean(V.next.temp)
	#V.last[,l]<-V.next.temp/mean(V.next.temp^2)
	}
	
	V.last[is.na(V.last)]<-0
	}


	U.next<-U.last
	V.next<-V.last
	
	#Uncommenting the next two lines eliminates the U and V
	#Gibbs update.
	#U.next<-svd0$u
	#V.next<-svd0$v

	Theta.last.0<-U.next%*%diag(D.post[1:min(n,k)])%*%t(V.next)
	Theta.last<-Theta.last.0+mean.mat
	#Not sure what the next three lines did.
	#chisq.ran<-rchisq(1,n*k)#+rchisq(min(n,k),1)
	#sigma<-(sum((Z.next-Theta.last)^2))/sum(chisq.ran)
	#Theta.last<-Theta.last
	Theta.mode<- U.next%*%diag(D.trunc[1:min(n,k)])%*%t(V.next)+mean.mat
	##Update muprime, invTau2, lambda.lasso
	muprime<-(abs(lambda.lasso*sqrt(sigma)/(abs.D.post)))
	invTau2<-sapply(1:min(n,k), FUN=function(i) rinv.gaussian(1, muprime[i], (lambda.lasso^2 ) ) )
	Dtau<-abs(1/invTau2)
	lambda.lasso<-rgamma(1, shape=min(n,k)+1  , rate=sum(Dtau[1:min(n,k)])/2+1.78  )^.5
	lambda.shrink<-(min(n,k)/(sum(Dtau[1:min(n,k)])/2+1.78))^.5
	
	
	return(list(
	"Theta.last"=Theta.last,
	"U.next"=U.next,
	"V.next"=V.next,
	"lambda.lasso"=lambda.lasso,
	"lambda.shrink"=lambda.shrink,
	"D.trunc"=D.trunc,
	"D.post"=D.post,
	"Theta.mode"=Theta.mode,
	"svd0"=svd0, 
	"Dtau"=Dtau,
	"D.ols"=svd0$d
	))
}



	expit<-function(x) exp(x)/(1+exp(x))
